# Install packages
# !pip install numpy
# !pip install matplotlib
# !pip install qiskit
# !pip install qiskit[visualization]
# !pip install qiskit[optimization]
# !pip install qiskit[machine-learning]
# Load packages
## General tools
import numpy as np
import matplotlib.pyplot as plt
from math import pi
## Qiskit Circuit Functions
from qiskit import *
from qiskit.quantum_info import *
from qiskit.visualization import *
We will be constructing Grover's Algorithm (also known as Grover's search algorithm).
For the examples in this textbook, our 'database' is comprised of all the possible computational basis states our qubits can be in. For example, if we have 3 qubits, our list is the states $|000\rangle, |001\rangle, \dots |111\rangle$ (i.e the states $|0\rangle \rightarrow |7\rangle$).
Oracle can be considered as an 'all knowing device" but we do not access to its informations or structures.
Grover’s algorithm solves oracles that add a negative phase to the solution states. I.e. for any state $|x\rangle$ in the computational basis:
If we initialize the 'output' qubit in the state $|{-}\rangle$, the phase kickback effect turns this into a Grover oracle (similar to the workings of the Deutsch-Jozsa oracle):
We then ignore the auxiliary ($|{-}\rangle$) qubit.
For the next part of this chapter, we aim to teach the core concepts of the algorithm. We will create example oracles where we know $\omega$ beforehand, and not worry ourselves with whether these oracles are useful or not.
So how does the algorithm work? Before looking at the list of items, we have no idea where the marked item is. Therefore, any guess of its location is as good as any other, which can be expressed in terms of a uniform superposition: $|s \rangle = \frac{1}{\sqrt{N}} \sum_{x = 0}^{N -1} | x\rangle.$
If at this point we were to measure in the standard basis $\{ | x \rangle \}$, this superposition would collapse, according to the fifth quantum law, to any one of the basis states with the same probability of $\frac{1}{N} = \frac{1}{2^n}$. Our chances of guessing the right value $w$ is therefore $1$ in $2^n$, as could be expected. Hence, on average we would need to try about $N/2 = 2^{n-1}$ times to guess the correct item.
Enter the procedure called amplitude amplification, which is how a quantum computer significantly enhances this probability. This procedure stretches out (amplifies) the amplitude of the marked item, which shrinks the other items' amplitude, so that measuring the final state will return the right item with near-certainty.
This algorithm has a nice geometrical interpretation in terms of two reflections, which generate a rotation in a two-dimensional plane. The only two special states we need to consider are the winner $| w \rangle$ and the uniform superposition $| s \rangle$. These two vectors span a two-dimensional plane in the vector space $\mathbb{C}^N.$ They are not quite perpendicular because $| w \rangle$ occurs in the superposition with amplitude $N^{-1/2}$ as well. We can, however, introduce an additional state $|s'\rangle$ that is in the span of these two vectors, which is perpendicular to $| w \rangle$ and is obtained from $|s \rangle$ by removing $| w \rangle$ and rescaling.
The left graphic corresponds to the two-dimensional plane spanned by perpendicular vectors $|w\rangle$ and $|s'\rangle$ which allows to express the initial state as $|s\rangle = \sin \theta | w \rangle + \cos \theta | s' \rangle,$ where $\theta = \arcsin \langle s | w \rangle = \arcsin \frac{1}{\sqrt{N}}$. The right graphic is a bar graph of the amplitudes of the state $| s \rangle$.
Geometrically this corresponds to a reflection of the state $|s\rangle$ about $|s'\rangle$. This transformation means that the amplitude in front of the $|w\rangle$ state becomes negative, which in turn means that the average amplitude (indicated by a dashed line) has been lowered.
Two reflections always correspond to a rotation. The transformation $U_s U_f$ rotates the initial state $|s\rangle$ closer towards the winner $|w\rangle$. The action of the reflection $U_s$ in the amplitude bar diagram can be understood as a reflection about the average amplitude. Since the average amplitude has been lowered by the first reflection, this transformation boosts the negative amplitude of $|w\rangle$ to roughly three times its original value, while it decreases the other amplitudes. We then go to step 2 to repeat the application. This procedure will be repeated several times to zero in on the winner.
After $t$ steps we will be in the state $|\psi_t\rangle$ where: $| \psi_t \rangle = (U_s U_f)^t | s \rangle.$
How many times do we need to apply the rotation? It turns out that roughly $\sqrt{N}$ rotations suffice. This becomes clear when looking at the amplitudes of the state $| \psi \rangle$. We can see that the amplitude of $| w \rangle$ grows linearly with the number of applications $\sim t N^{-1/2}$. However, since we are dealing with amplitudes and not probabilities, the vector space's dimension enters as a square root. Therefore it is the amplitude, and not just the probability, which is being amplified in this procedure.
In the case that there are multiple solutions, $M$, it can be shown that roughly $\sqrt{(N/M)}$ rotations will suffice.
Let's first have a look at the case of Grover's algorithm for $N=4$ which is realized with 2 qubits. In this particular case, only one rotation is required to rotate the initial state $|s\rangle$ to the winner $|w\rangle$:
1) Following the above introduction, in the case $N=4$ we have $$\theta = \arcsin \frac{1}{2} = \frac{\pi}{6}.$$
2) After $t$ steps, we have $$(U_s U_\omega)^t | s \rangle = \sin \theta_t | \omega \rangle + \cos \theta_t | s' \rangle ,$$where $$\theta_t = (2t+1)\theta.$$
3) In order to obtain $| \omega \rangle$ we need $\theta_t = \frac{\pi}{2}$, which with $\theta=\frac{\pi}{6}$ inserted above results to $t=1$. This implies that after $t=1$ rotation the searched element is found.
We will now follow through an example using a specific oracle.
Oracle for $\lvert \omega \rangle = \lvert 11 \rangle$
or:
$$ U_\omega = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & -1 \\ \end{bmatrix} $$which you may recognize as the controlled-Z gate. Thus, for this example, our oracle is simply the controlled-Z gate:
Reflection $U_s$
In order to complete the circuit we need to implement the additional reflection $U_s = 2|s\rangle\langle s| - \mathbb{1}$. Since this is a reflection about $|s\rangle$, we want to add a negative phase to every state orthogonal to $|s\rangle$.
One way we can do this is to use the operation that transforms the state $|s\rangle \rightarrow |0\rangle$, which we already know is the Hadamard gate applied to each qubit:
i.e. the signs of each state are flipped except for $\lvert 00 \rangle$. As can easily be verified, one way of implementing $U_0$ is the following circuit:
The complete circuit for $U_s$ looks like this:
Full Circuit for $\lvert w \rangle = |11\rangle$ Since in the particular case of $N=4$ only one rotation is required we can combine the above components to build the full circuit for Grover's algorithm for the case $\lvert w \rangle = |11\rangle$:
# Initialization
import matplotlib.pyplot as plt
import numpy as np
# Importing Qiskit
from qiskit import IBMQ, Aer, assemble, transpile
from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister
from qiskit.providers.ibmq import least_busy
# Import basic plot tools
from qiskit.visualization import plot_histogram
n = 2
grover_circuit = QuantumCircuit(n)
def initialize_s(qc, qubits):
"""Apply a H-gate to 'qubits' in qc"""
for q in qubits:
qc.h(q)
return qc
grover_circuit = initialize_s(grover_circuit, [0,1])
grover_circuit.barrier()
grover_circuit.draw(output='mpl')
grover_circuit.cz(0,1) # Oracle
grover_circuit.barrier()
grover_circuit.draw(output='mpl')
# Diffusion operator (U_s)
grover_circuit.h([0,1])
grover_circuit.z([0,1])
grover_circuit.cz(0,1)
grover_circuit.h([0,1])
grover_circuit.barrier()
grover_circuit.draw(output='mpl')
Let's run the circuit in simulation. First, we can verify that we have the correct statevector:
sim = Aer.get_backend('aer_simulator')
# we need to make a copy of the circuit with the 'save_statevector'
# instruction to run on the Aer simulator
grover_circuit_sim = grover_circuit.copy()
grover_circuit_sim.save_statevector()
qobj = assemble(grover_circuit_sim)
result = sim.run(qobj).result()
statevec = result.get_statevector()
# Display the output state as a vector
display(array_to_latex(statevec, prefix="\\text{Statevector } |\\psi\\rangle = "))
# from qiskit_textbook.tools import vector2latex
# vector2latex(statevec, pretext="|\\psi\\rangle =")
grover_circuit.measure_all()
aer_sim = Aer.get_backend('aer_simulator')
qobj = assemble(grover_circuit)
result = aer_sim.run(qobj).result()
counts = result.get_counts()
plot_histogram(counts)
We now go through the example of Grover's algorithm for 3 qubits with two marked states $\lvert101\rangle$ and $\lvert110\rangle$, following the implementation found in Reference [2]. The quantum circuit to solve the problem using a phase oracle is:
1) Apply Hadamard gates to $3$ qubits initialized to $\lvert000\rangle$ to create a uniform superposition: $$\lvert \psi_1 \rangle = \frac{1}{\sqrt{8}} \left( \lvert000\rangle + \lvert001\rangle + \lvert010\rangle + \lvert011\rangle + \lvert100\rangle + \lvert101\rangle + \lvert110\rangle + \lvert111\rangle \right) $$
2) Mark states $\lvert101\rangle$ and $\lvert110\rangle$ using a phase oracle: $$\lvert \psi_2 \rangle = \frac{1}{\sqrt{8}} \left( \lvert000\rangle + \lvert001\rangle + \lvert010\rangle + \lvert011\rangle + \lvert100\rangle - \lvert101\rangle - \lvert110\rangle + \lvert111\rangle \right) $$
3) Perform the reflection around the average amplitude:
A. Apply Hadamard gates to the qubits $$\lvert \psi_{3a} \rangle = \frac{1}{2} \left( \lvert000\rangle +\lvert011\rangle +\lvert100\rangle -\lvert111\rangle \right) $$
B. Apply X gates to the qubits $$\lvert \psi_{3b} \rangle = \frac{1}{2} \left( -\lvert000\rangle +\lvert011\rangle +\lvert100\rangle +\lvert111\rangle \right) $$
C. Apply a doubly controlled Z gate between the 1, 2 (controls) and 3 (target) qubits $$\lvert \psi_{3c} \rangle = \frac{1}{2} \left( -\lvert000\rangle +\lvert011\rangle +\lvert100\rangle -\lvert111\rangle \right) $$
D. Apply X gates to the qubits $$\lvert \psi_{3d} \rangle = \frac{1}{2} \left( -\lvert000\rangle +\lvert011\rangle +\lvert100\rangle -\lvert111\rangle \right) $$
E. Apply Hadamard gates to the qubits $$\lvert \psi_{3e} \rangle = \frac{1}{\sqrt{2}} \left( -\lvert101\rangle -\lvert110\rangle \right) $$
4) Measure the $3$ qubits to retrieve states $\lvert101\rangle$ and $\lvert110\rangle$
Note that since there are 2 solutions and 8 possibilities, we will only need to run one iteration (steps 2 & 3).
qc = QuantumCircuit(3)
qc.cz(0, 1)
oracle_ex3 = qc.to_gate()
oracle_ex3.name = "U$_\omega$"
Details: Creating a General Diffuser
Remember that we can create $U_s$ from $U_0$:
$$ U_s = H^{\otimes n} U_0 H^{\otimes n} $$And a multi-controlled-Z gate ($MCZ$) inverts the phase of the state $|11\dots 1\rangle$:
$$ MCZ = \begin{bmatrix} 1 & 0 & 0 & \cdots & 0 \\ 0 & 1 & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & -1 \\ \end{bmatrix} \begin{aligned} \\ \\ \\ \leftarrow \text{Add negative phase to} \; |11\dots 1\rangle\\ \end{aligned} $$Applying an X-gate to each qubit performs the transformation:
$$ \begin{aligned} |00\dots 0\rangle & \rightarrow |11\dots 1\rangle\\ |11\dots 1\rangle & \rightarrow |00\dots 0\rangle \end{aligned} $$So:
$$ U_0 = - X^{\otimes n} (MCZ) X^{\otimes n} $$Using these properties together, we can create $U_s$ using H-gates, X-gates, and a single multi-controlled-Z gate:
$$ U_s = - H^{\otimes n} U_0 H^{\otimes n} = H^{\otimes n} X^{\otimes n} (MCZ) X^{\otimes n} H^{\otimes n} $$Note that we can ignore the global phase of -1.
def diffuser(nqubits):
qc = QuantumCircuit(nqubits)
# Apply transformation |s> -> |00..0> (H-gates)
for qubit in range(nqubits):
qc.h(qubit)
# Apply transformation |00..0> -> |11..1> (X-gates)
for qubit in range(nqubits):
qc.x(qubit)
# Do multi-controlled-Z gate
qc.h(nqubits-1)
qc.mct(list(range(nqubits-1)), nqubits-1) # multi-controlled-toffoli
qc.h(nqubits-1)
# Apply transformation |11..1> -> |00..0>
for qubit in range(nqubits):
qc.x(qubit)
# Apply transformation |00..0> -> |s>
for qubit in range(nqubits):
qc.h(qubit)
# We will return the diffuser as a gate
U_s = qc.to_gate()
U_s.name = "U$_s$"
return U_s
n = 3
grover_circuit = QuantumCircuit(n)
grover_circuit = initialize_s(grover_circuit, [0,1,2])
grover_circuit.barrier()
grover_circuit.append(oracle_ex3, [0,1,2])
grover_circuit.barrier()
grover_circuit.append(diffuser(n), [0,1,2])
grover_circuit.measure_all()
grover_circuit.draw(output='mpl')
decompose() on circuit object.grover_circuit.decompose().draw(output='mpl')
aer_sim = Aer.get_backend('aer_simulator')
transpiled_grover_circuit = transpile(grover_circuit, aer_sim)
qobj = assemble(transpiled_grover_circuit)
results = aer_sim.run(qobj).result()
counts = results.get_counts()
plot_histogram(counts)
Regression using quantum codes:
EstimatorQNNVQR)Link: Neural Network Regressor
IBMQ.load_account()
<AccountProvider for IBMQ(hub='ibm-q', group='open', project='main')>
from qiskit import QuantumCircuit
from qiskit.algorithms.optimizers import COBYLA, L_BFGS_B
from qiskit.circuit import Parameter
from qiskit.circuit.library import RealAmplitudes, ZZFeatureMap
from qiskit.utils import algorithm_globals
from qiskit_machine_learning.algorithms.classifiers import NeuralNetworkClassifier, VQC
from qiskit_machine_learning.algorithms.regressors import NeuralNetworkRegressor, VQR
from qiskit_machine_learning.neural_networks import SamplerQNN, EstimatorQNN
from qiskit.utils import QuantumInstance, algorithm_globals
from qiskit.circuit import Parameter
# from qiskit.circuit.library import RealAmplitudes, ZZFeatureMap
from qiskit.algorithms.optimizers import COBYLA, L_BFGS_B
from qiskit_machine_learning.neural_networks import TwoLayerQNN, CircuitQNN
# from qiskit_machine_learning.algorithms.classifiers import NeuralNetworkClassifier, VQC
from qiskit_machine_learning.algorithms.regressors import NeuralNetworkRegressor, VQR
from IPython.display import clear_output
# Loading Quantum backend - aer_simulator
quantum_instance = QuantumInstance(Aer.get_backend("aer_simulator"), shots=1024)
algorithm_globals.random_seed = 42 # setting the random_seed number
num_samples = 20 # the number of observations of the data to create
eps = 0.2 # error coefficient to be multiplicated
lb, ub = -np.pi, np.pi # lower bound = -pi, upper bound = +pi
# Create a dataset generated from the sine function (true function values)
X_ = np.linspace(lb, ub, num=50).reshape(50, 1)
f = lambda x: np.sin(x)
# Adding noise values to x and y values from the true function values (sine function values)
X = (ub - lb) * algorithm_globals.random.random([num_samples, 1]) + lb
y = f(X[:, 0]) + eps * (2 * algorithm_globals.random.random(num_samples) - 1)
# Datasets
X
array([[ 1.72131662],
[-0.38403809],
[ 2.25313718],
[ 1.2400999 ],
[-2.54985893],
[ 2.98842337],
[ 1.64078914],
[ 1.79739504],
[-2.33663096],
[-0.31173435],
[-0.81179996],
[ 2.68144351],
[ 0.90393121],
[ 2.02797103],
[-0.35553907],
[-1.71380966],
[ 0.34296633],
[-2.74061701],
[ 2.05856737],
[ 0.82727182]])
y
array([ 1.09192829, -0.43285708, 0.96437789, 1.10306489, -0.44644741,
0.03042655, 0.9842399 , 0.79195768, -0.8590883 , -0.23349027,
-0.6276222 , 0.63108564, 0.71609467, 0.84548709, -0.36027345,
-1.11400247, 0.18825075, -0.40003482, 0.77414341, 0.80401304])
# Plotting true function values (sine function) and the generated data
plt.figure(figsize=(8,6))
plt.plot(X_, f(X_), "r--")
plt.plot(X, y, "bo")
plt.show()
Red dotted line: True function value $f(x) = sin(x)$, where $x \in (-\pi, \pi)$Blue dots: generated data from true function values added with noisy valuesred dotted line) with only having / using the data we have, which is blue dots.EstimatorQNN¶EstimatorQNN that returns values in $\big[-1, 1\big]$.SamplerQNN but that exceeds the scope of this tutorial.# construct simple feature map
param_x = Parameter("x")
feature_map = QuantumCircuit(1, name="fm")
feature_map.ry(param_x, 0)
feature_map.draw(output = 'mpl')
# construct simple ansatz
param_y = Parameter("y")
ansatz = QuantumCircuit(1, name="vf")
ansatz.ry(param_y, 0)
ansatz.draw(output = 'mpl')
# construct a circuit
qc = QuantumCircuit(1)
qc.compose(feature_map, inplace=True)
qc.compose(ansatz, inplace=True)
qc.draw(output = 'mpl')
# construct QNN
regression_estimator_qnn = EstimatorQNN(
circuit=qc, input_params=feature_map.parameters, weight_params=ansatz.parameters
)
# callback function that draws a live plot when the .fit() method is called
def callback_graph(weights, obj_func_eval):
clear_output(wait=True)
objective_func_vals.append(obj_func_eval)
plt.title("Objective function value against iteration")
plt.xlabel("Iteration")
plt.ylabel("Objective function value")
plt.plot(range(len(objective_func_vals)), objective_func_vals)
plt.show()
# construct the regressor from the neural network
# callback function that draws a live plot when the .fit() method is called
regressor = NeuralNetworkRegressor(
neural_network=regression_estimator_qnn,
loss="squared_error",
optimizer=L_BFGS_B(maxiter=5),
callback=callback_graph,
)
# create empty array for callback to store evaluations of the objective function
objective_func_vals = []
plt.rcParams["figure.figsize"] = (12, 6)
# fit to data
regressor.fit(X, y)
# return to default figsize
plt.rcParams["figure.figsize"] = (6, 4)
# score the result - returns the mean accuracy of the score (R^2)
regressor.score(X, y)
0.9733353285953001
# plot target function
plt.figure(figsize=(8,6))
plt.plot(X_, f(X_), "r--")
# plot data
plt.plot(X, y, "bo")
# plot fitted line
y_ = regressor.predict(X_)
plt.plot(X_, y_, "g-")
plt.show()
VQR is a special variant of the NeuralNetworkRegressor with a OpflowQNN.
By default it considers the L2Loss function to minimize the mean squared error between predictions and targets.
vqr = VQR(
feature_map=feature_map,
ansatz=ansatz,
optimizer=L_BFGS_B(),
quantum_instance=quantum_instance,
callback=callback_graph,
)
C:\Users\user\AppData\Local\Temp\ipykernel_9996\1717641275.py:1: DeprecationWarning: The quantum_instance argument is deprecated as of version 0.5.0 and will be removed no sooner than 3 months after the release. Instead use the estimator argument. vqr = VQR(
# create empty array for callback to store evaluations of the objective function
objective_func_vals = []
plt.rcParams["figure.figsize"] = (12, 6)
# fit regressor
vqr.fit(X, y)
# return to default figsize
plt.rcParams["figure.figsize"] = (6, 4)
# score the result - returns the mean accuracy of the score (R^2)
vqr.score(X, y)
0.9690894881891703
# plot target function
plt.figure(figsize=(8,6))
plt.plot(X_, f(X_), "r--")
# plot data
plt.plot(X, y, "bo")
# plot fitted line
y_ = vqr.predict(X_)
plt.plot(X_, y_, "g-")
plt.show()
#
Various types of quantum computers
For more information, here are some useful links to check.
To install PennyLane
Link: PennyLane installation
# install the latest released version of PennyLane
#!pip install pennylane
from pennylane import numpy as np
import matplotlib.pyplot as plt
from math import pi
import pennylane as qml
We will learn how to build quantum circuits using PennyLane codes.
Defining a device
To run—and later optimize—a quantum circuit, one needs to first specify a computational device.
device(name, wires=1, *args, **kwargs)
name: the name of the device / simulators to load.wires: the number of wires to initialize. Think the wire as the number of qubits.shots: the number of execution of the circuit.# Setup the device using qml.device()
## Assigning the device and defining the measurement of the circuit
dev = qml.device(name = 'default.qubit', wires = 2)
@qml.qnode(dev)
def circuit():
return qml.expval(qml.PauliZ(0) @ qml.PauliZ(1))
drawer = qml.draw(circuit, show_all_wires=True)
print(drawer())
0: ───┤ ╭<Z@Z> 1: ───┤ ╰<Z@Z>
# Setup the device using qml.device()
## Assigning names to each wires
dev = qml.device(name = 'default.qubit', wires = ['a', 'b'])
@qml.qnode(dev)
def circuit(x):
return qml.expval(qml.PauliZ('a') @ qml.PauliZ('b'))
drawer = qml.draw(circuit, show_all_wires=True)
print(drawer(0.5))
a: ───┤ ╭<Z@Z> b: ───┤ ╰<Z@Z>
# Setup the device using qml.device()
## Assigning the number of shots
shots_list = [5, 10, 1000] # number of shots for each batches
dev = qml.device("default.qubit", wires=2, shots=shots_list)
@qml.qnode(dev)
def circuit(x):
qml.RX(x, wires=0)
qml.CNOT(wires=[0, 1])
return qml.expval(qml.PauliZ(0) @ qml.PauliX(1)), qml.expval(qml.PauliZ(0))
drawer = qml.draw(circuit, show_all_wires=True)
print(drawer(x = 0.5))
0: ──RX(0.50)─╭●─┤ ╭<Z@X> <Z> 1: ───────────╰X─┤ ╰<Z@X>
## The result of the execution of the circuit.
circuit(0.5)
tensor([[-0.2 , 1. ],
[-0.6 , 1. ],
[ 0.002, 0.876]], requires_grad=True)
Creating a quantum node (Quantum operators + measurement)
dev = qml.device('default.qubit', wires=['aux', 'q1', 'q2'])
# Quntum function
def my_quantum_function(x, y):
qml.RZ(x, wires='q1')
qml.CNOT(wires=['aux' ,'q1'])
qml.RY(y, wires='q2')
return qml.expval(qml.PauliZ('q2'))
circuit = qml.QNode(my_quantum_function, dev)
circuit(x = np.pi/4, y = 0.7)
tensor(0.76484219, requires_grad=True)
# Plot the circuit
print(qml.draw(circuit)(np.pi/4, 0.7))
aux: ───────────╭●─┤ q1: ──RZ(0.79)─╰X─┤ q2: ──RY(0.70)────┤ <Z>
# Plot the circuit with matplotlib.pyplot library
import matplotlib.pyplot as plt
fig, ax = qml.draw_mpl(circuit)(x = np.pi/4, y = 0.7)
plt.show()
Measurement
# Measurement 1) Sampling result from each shots
dev = qml.device("default.qubit", wires=2, shots=1000)
@qml.qnode(dev)
def circuit():
qml.Hadamard(wires=0)
qml.CNOT(wires=[0, 1])
return qml.sample(qml.PauliZ(0)), qml.sample(qml.PauliZ(1))
fig, ax = qml.draw_mpl(circuit)()
plt.show()
result_sample = circuit()
# Print out the dimension of the result
print(result_sample.shape)
# Let's look at the result
result_sample
(2, 1000)
tensor([[-1, -1, 1, ..., 1, 1, 1],
[-1, -1, 1, ..., 1, 1, 1]], dtype=int64, requires_grad=True)
1 corresponds to $|0\rangle$, and -1 corresponds to $|1\rangle$.# Measurement 2) Sampling result from each shots
dev = qml.device("default.qubit", wires=2, shots=1000)
@qml.qnode(dev)
def circuit():
qml.Hadamard(wires=0)
qml.CNOT(wires=[0, 1])
return qml.probs(wires=[0, 1])
fig, ax = qml.draw_mpl(circuit)()
plt.show()
result_prob = circuit()
# Print out the dimension of the result
print(result_prob.shape)
# Let's look at the result
result_prob
(4,)
tensor([0.471, 0. , 0. , 0.529], requires_grad=True)
# Measurement 3) Sampling result from each shots
dev = qml.device("default.qubit", wires=3, shots=1000)
@qml.qnode(dev)
def my_quantum_function(x, y):
qml.RZ(x, wires=0)
qml.CNOT(wires=[0, 1])
qml.RY(y, wires=1)
qml.CNOT(wires=[0, 2])
return qml.expval(qml.PauliZ(0) @ qml.Identity(1) @qml.PauliX(2))
fig, ax = qml.draw_mpl(my_quantum_function)(x = pi/2, y = pi/3)
plt.show()
result_tensor = my_quantum_function(x = pi/2, y = pi/3)
# Let's look at the result
result_tensor
tensor(-0.004, requires_grad=True)
This optimization example demonstrates how to encode binary inputs into the initial state of the variational circuit, which is simply a computational basis state.
We then show how to encode real vectors as amplitude vectors (amplitude encoding) and train the model to recognize the first two classes of flowers in the Iris dataset.
Load packages / functions
import pennylane as qml
from pennylane import numpy as np
from pennylane.optimize import NesterovMomentumOptimizer
dev = qml.device("default.qubit", wires=4)
Variational classifiers usually define a “layer” or “block”, which is an elementary circuit architecture that gets repeated to build the variational circuit.
Our circuit layer consists of an arbitrary rotation on every qubit, as well as CNOTs that entangle each qubit with its neighbor.
qml.Rot is a code function that performs an arbitrary rotational gate operator.
def layer(W):
qml.Rot(W[0, 0], W[0, 1], W[0, 2], wires=0)
qml.Rot(W[1, 0], W[1, 1], W[1, 2], wires=1)
qml.Rot(W[2, 0], W[2, 1], W[2, 2], wires=2)
qml.Rot(W[3, 0], W[3, 1], W[3, 2], wires=3)
qml.CNOT(wires=[0, 1])
qml.CNOT(wires=[1, 2])
qml.CNOT(wires=[2, 3])
qml.CNOT(wires=[3, 0])
We also need a way to encode data inputs $x$ into the circuit, so that the measured output depends on the inputs. In this first example, the inputs are bitstrings, which we encode into the state of the qubits. The quantum state $\psi$ after state preparation is a computational basis state that has 1s where $x$ has 1s, for example,
$$ x = 0101 \rightarrow \vert\psi\rangle = \vert0101\rangle. $$We use the BasisState function provided by PennyLane, which expects $x$ to be a list of zeros and ones, i.e. [0,1,0,1].
def statepreparation(x):
qml.BasisState(x, wires=[0, 1, 2, 3])
Now we define the quantum node as a state preparation routine, followed by a repetition of the layer structure. Borrowing from machine learning, we call the parameters weights.
@qml.qnode(dev)
def circuit(weights, x):
statepreparation(x)
for W in weights:
layer(W)
return qml.expval(qml.PauliZ(0))
Different from previous examples, the quantum node takes the data as a keyword argument $x$ (with the default value None). Keyword arguments of a quantum node are considered as fixed when calculating a gradient; they are never trained.
If we want to add a “classical” bias parameter, the variational quantum classifier also needs some post-processing. We define the final model by a classical node that uses the first variable, and feeds the remainder into the quantum node. Before this, we reshape the list of remaining variables for easy use in the quantum node.
def variational_classifier(weights, bias, x):
return circuit(weights, x) + bias
def square_loss(labels, predictions):
loss = 0
for l, p in zip(labels, predictions):
loss = loss + (l - p) ** 2
loss = loss / len(labels)
return loss
def accuracy(labels, predictions):
loss = 0
for l, p in zip(labels, predictions):
if abs(l - p) < 1e-5:
loss = loss + 1
loss = loss / len(labels)
return loss
def cost(weights, bias, X, Y):
predictions = [variational_classifier(weights, bias, x) for x in X]
return square_loss(Y, predictions)
from pennylane import numpy as np
data = np.array([[0., 0., 0., 0., 0],
[0., 0., 0., 1., 1],
[0., 0., 1., 0., 1],
[0., 0., 1., 1., 0],
[0., 1., 0., 0., 1],
[0., 1., 0., 1., 0],
[0., 1., 1., 0., 0],
[0., 1., 1., 1., 1],
[1., 0., 0., 0., 1],
[1., 0., 0., 1., 0],
[1., 0., 1., 0., 0],
[1., 0., 1., 1., 1],
[1., 1., 0., 0., 0],
[1., 1., 0., 1., 1],
[1., 1., 1., 0., 1],
[1., 1., 1., 1., 0]])
X = np.array(data[:, :-1], requires_grad=False)
Y = np.array(data[:, -1], requires_grad=False)
Y = Y * 2 - np.ones(len(Y)) # shift label from {0, 1} to {-1, 1}
for i in range(5):
print("X = {}, Y = {: d}".format(X[i], int(Y[i])))
print("...")
X = [0. 0. 0. 0.], Y = -1 X = [0. 0. 0. 1.], Y = 1 X = [0. 0. 1. 0.], Y = 1 X = [0. 0. 1. 1.], Y = -1 X = [0. 1. 0. 0.], Y = 1 ...
np.random.seed(0)
num_qubits = 4
num_layers = 2
weights_init = 0.01 * np.random.randn(num_layers, num_qubits, 3, requires_grad=True)
bias_init = np.array(0.0, requires_grad=True)
print(weights_init, bias_init)
[[[ 0.01764052 0.00400157 0.00978738] [ 0.02240893 0.01867558 -0.00977278] [ 0.00950088 -0.00151357 -0.00103219] [ 0.00410599 0.00144044 0.01454274]] [[ 0.00761038 0.00121675 0.00443863] [ 0.00333674 0.01494079 -0.00205158] [ 0.00313068 -0.00854096 -0.0255299 ] [ 0.00653619 0.00864436 -0.00742165]]] 0.0
Next we create an optimizer and choose a batch size…
opt = NesterovMomentumOptimizer(0.5)
batch_size = 5
…and train the optimizer. We track the accuracy - the share of correctly classified data samples. For this we compute the outputs of the variational classifier and turn them into predictions in $\big\{-1, 1\big\}$ by taking the sign of the output.
weights = weights_init
bias = bias_init
iter_num = 10
for it in range(iter_num):
# Update the weights by one optimizer step
batch_index = np.random.randint(0, len(X), (batch_size,))
X_batch = X[batch_index]
Y_batch = Y[batch_index]
weights, bias, _, _ = opt.step(cost, weights, bias, X_batch, Y_batch)
# Compute accuracy
predictions = [np.sign(variational_classifier(weights, bias, x)) for x in X]
acc = accuracy(Y, predictions)
print(
"Iter: {:5d} | Cost: {:0.7f} | Accuracy: {:0.7f} ".format(
it + 1, cost(weights, bias, X, Y), acc
)
)
Iter: 1 | Cost: 3.4355534 | Accuracy: 0.5000000 Iter: 2 | Cost: 1.9717733 | Accuracy: 0.5000000 Iter: 3 | Cost: 1.8182812 | Accuracy: 0.5000000 Iter: 4 | Cost: 1.5042404 | Accuracy: 0.5000000 Iter: 5 | Cost: 1.1477739 | Accuracy: 0.5000000 Iter: 6 | Cost: 1.2734990 | Accuracy: 0.6250000 Iter: 7 | Cost: 0.8290628 | Accuracy: 0.5000000 Iter: 8 | Cost: 0.3226183 | Accuracy: 1.0000000 Iter: 9 | Cost: 0.1436206 | Accuracy: 1.0000000 Iter: 10 | Cost: 0.2982810 | Accuracy: 1.0000000
To encode real-valued vectors into the amplitudes of a quantum state, we use a 2-qubit simulator.
dev = qml.device("default.qubit", wires=2)
def get_angles(x):
beta0 = 2 * np.arcsin(np.sqrt(x[1] ** 2) / np.sqrt(x[0] ** 2 + x[1] ** 2 + 1e-12))
beta1 = 2 * np.arcsin(np.sqrt(x[3] ** 2) / np.sqrt(x[2] ** 2 + x[3] ** 2 + 1e-12))
beta2 = 2 * np.arcsin(
np.sqrt(x[2] ** 2 + x[3] ** 2)
/ np.sqrt(x[0] ** 2 + x[1] ** 2 + x[2] ** 2 + x[3] ** 2)
)
return np.array([beta2, -beta1 / 2, beta1 / 2, -beta0 / 2, beta0 / 2])
def statepreparation(a):
qml.RY(a[0], wires=0)
qml.CNOT(wires=[0, 1])
qml.RY(a[1], wires=1)
qml.CNOT(wires=[0, 1])
qml.RY(a[2], wires=1)
qml.PauliX(wires=0)
qml.CNOT(wires=[0, 1])
qml.RY(a[3], wires=1)
qml.CNOT(wires=[0, 1])
qml.RY(a[4], wires=1)
qml.PauliX(wires=0)
x = np.array([0.53896774, 0.79503606, 0.27826503, 0.0], requires_grad=False)
ang = get_angles(x)
@qml.qnode(dev)
def test(angles):
statepreparation(angles)
return qml.expval(qml.PauliZ(0))
test(ang)
print("x : ", x)
print("angles : ", ang)
print("amplitude vector: ", np.real(dev.state))
x : [0.53896774 0.79503606 0.27826503 0. ] angles : [ 0.56397465 -0. 0. -0.97504604 0.97504604] amplitude vector: [ 5.38967743e-01 7.95036065e-01 2.78265032e-01 -2.77555756e-17]
Note that the default.qubit simulator provides a shortcut to statepreparation with the command qml.QubitStateVector(x, wires=[0, 1]). However, some devices may not support an arbitrary state-preparation routine.
Since we are working with only 2 qubits now, we need to update the layer function as well.
def layer(W):
qml.Rot(W[0, 0], W[0, 1], W[0, 2], wires=0)
qml.Rot(W[1, 0], W[1, 1], W[1, 2], wires=1)
qml.CNOT(wires=[0, 1])
The variational classifier model and its cost remain essentially the same, but we have to reload them with the new state preparation and layer functions.
@qml.qnode(dev)
def circuit(weights, angles):
statepreparation(angles)
for W in weights:
layer(W)
return qml.expval(qml.PauliZ(0))
def variational_classifier(weights, bias, angles):
return circuit(weights, angles) + bias
def cost(weights, bias, features, labels):
predictions = [variational_classifier(weights, bias, f) for f in features]
return square_loss(labels, predictions)
# !pip install pandas
# !pip install sklearn
from sklearn.datasets import load_iris
import pandas as pd
import numpy as np
# visualization package
import matplotlib.pyplot as plt
import seaborn as sns
iris = load_iris() # sample data load
#print(iris) # pring out data with variable types and its description
print(iris.DESCR) # Description of the dataset
.. _iris_dataset:
Iris plants dataset
--------------------
**Data Set Characteristics:**
:Number of Instances: 150 (50 in each of three classes)
:Number of Attributes: 4 numeric, predictive attributes and the class
:Attribute Information:
- sepal length in cm
- sepal width in cm
- petal length in cm
- petal width in cm
- class:
- Iris-Setosa
- Iris-Versicolour
- Iris-Virginica
:Summary Statistics:
============== ==== ==== ======= ===== ====================
Min Max Mean SD Class Correlation
============== ==== ==== ======= ===== ====================
sepal length: 4.3 7.9 5.84 0.83 0.7826
sepal width: 2.0 4.4 3.05 0.43 -0.4194
petal length: 1.0 6.9 3.76 1.76 0.9490 (high!)
petal width: 0.1 2.5 1.20 0.76 0.9565 (high!)
============== ==== ==== ======= ===== ====================
:Missing Attribute Values: None
:Class Distribution: 33.3% for each of 3 classes.
:Creator: R.A. Fisher
:Donor: Michael Marshall (MARSHALL%PLU@io.arc.nasa.gov)
:Date: July, 1988
The famous Iris database, first used by Sir R.A. Fisher. The dataset is taken
from Fisher's paper. Note that it's the same as in R, but not as in the UCI
Machine Learning Repository, which has two wrong data points.
This is perhaps the best known database to be found in the
pattern recognition literature. Fisher's paper is a classic in the field and
is referenced frequently to this day. (See Duda & Hart, for example.) The
data set contains 3 classes of 50 instances each, where each class refers to a
type of iris plant. One class is linearly separable from the other 2; the
latter are NOT linearly separable from each other.
.. topic:: References
- Fisher, R.A. "The use of multiple measurements in taxonomic problems"
Annual Eugenics, 7, Part II, 179-188 (1936); also in "Contributions to
Mathematical Statistics" (John Wiley, NY, 1950).
- Duda, R.O., & Hart, P.E. (1973) Pattern Classification and Scene Analysis.
(Q327.D83) John Wiley & Sons. ISBN 0-471-22361-1. See page 218.
- Dasarathy, B.V. (1980) "Nosing Around the Neighborhood: A New System
Structure and Classification Rule for Recognition in Partially Exposed
Environments". IEEE Transactions on Pattern Analysis and Machine
Intelligence, Vol. PAMI-2, No. 1, 67-71.
- Gates, G.W. (1972) "The Reduced Nearest Neighbor Rule". IEEE Transactions
on Information Theory, May 1972, 431-433.
- See also: 1988 MLC Proceedings, 54-64. Cheeseman et al"s AUTOCLASS II
conceptual clustering system finds 3 classes in the data.
- Many, many more ...
# feature_names(x variables) 와 target(y variable)을 잘 나타내도록 data frame 생성
df = pd.DataFrame(data=iris.data, columns=iris.feature_names)
df['target'] = iris.target
# 숫자형으로 기록된 target(y variable) - (0.0, 1.0, 2.0)을 문자형으로 변환
df['target'] = df['target'].map({0:"setosa", 1:"versicolor", 2:"virginica"})
print(df)
sns.pairplot(df, hue="target", height=3)
plt.show()
sepal length (cm) sepal width (cm) petal length (cm) petal width (cm) \
0 5.1 3.5 1.4 0.2
1 4.9 3.0 1.4 0.2
2 4.7 3.2 1.3 0.2
3 4.6 3.1 1.5 0.2
4 5.0 3.6 1.4 0.2
.. ... ... ... ...
145 6.7 3.0 5.2 2.3
146 6.3 2.5 5.0 1.9
147 6.5 3.0 5.2 2.0
148 6.2 3.4 5.4 2.3
149 5.9 3.0 5.1 1.8
target
0 setosa
1 setosa
2 setosa
3 setosa
4 setosa
.. ...
145 virginica
146 virginica
147 virginica
148 virginica
149 virginica
[150 rows x 5 columns]
from pennylane import numpy as np
data = np.loadtxt("iris_classes1and2_scaled.txt")
X = data[:, 0:2]
print("First X sample (original) :", X[0])
# pad the vectors to size 2^2 with constant values
padding = 0.3 * np.ones((len(X), 1))
X_pad = np.c_[np.c_[X, padding], np.zeros((len(X), 1))]
print("First X sample (padded) :", X_pad[0])
# normalize each input
normalization = np.sqrt(np.sum(X_pad ** 2, -1))
X_norm = (X_pad.T / normalization).T
print("First X sample (normalized):", X_norm[0])
# angles for state preparation are new features
features = np.array([get_angles(x) for x in X_norm], requires_grad=False)
print("First features sample :", features[0])
Y = data[:, -1]
First X sample (original) : [0.4 0.75] First X sample (padded) : [0.4 0.75 0.3 0. ] First X sample (normalized): [0.44376016 0.83205029 0.33282012 0. ] First features sample : [ 0.67858523 -0. 0. -1.080839 1.080839 ]
import matplotlib.pyplot as plt
plt.figure()
plt.scatter(X[:, 0][Y == 1], X[:, 1][Y == 1], c="b", marker="o", edgecolors="k")
plt.scatter(X[:, 0][Y == -1], X[:, 1][Y == -1], c="r", marker="o", edgecolors="k")
plt.title("Original data")
plt.show()
plt.figure()
dim1 = 0
dim2 = 1
plt.scatter(
X_norm[:, dim1][Y == 1], X_norm[:, dim2][Y == 1], c="b", marker="o", edgecolors="k"
)
plt.scatter(
X_norm[:, dim1][Y == -1], X_norm[:, dim2][Y == -1], c="r", marker="o", edgecolors="k"
)
plt.title("Padded and normalised data (dims {} and {})".format(dim1, dim2))
plt.show()
plt.figure()
dim1 = 0
dim2 = 3
plt.scatter(
features[:, dim1][Y == 1], features[:, dim2][Y == 1], c="b", marker="o", edgecolors="k"
)
plt.scatter(
features[:, dim1][Y == -1], features[:, dim2][Y == -1], c="r", marker="o", edgecolors="k"
)
plt.title("Feature vectors (dims {} and {})".format(dim1, dim2))
plt.show()
np.random.seed(0)
num_data = len(Y)
num_train = int(0.75 * num_data)
index = np.random.permutation(range(num_data))
feats_train = features[index[:num_train]]
Y_train = Y[index[:num_train]]
feats_val = features[index[num_train:]]
Y_val = Y[index[num_train:]]
# We need these later for plotting
X_train = X[index[:num_train]]
X_val = X[index[num_train:]]
num_qubits = 2
num_layers = 6
weights_init = 0.01 * np.random.randn(num_layers, num_qubits, 3, requires_grad=True)
bias_init = np.array(0.0, requires_grad=True)
opt = NesterovMomentumOptimizer(0.01)
batch_size = 50
# train the variational classifier
weights = weights_init
bias = bias_init
Iter_num = 25
for it in range(Iter_num):
# Update the weights by one optimizer step
batch_index = np.random.randint(0, num_train, (batch_size,))
feats_train_batch = feats_train[batch_index]
Y_train_batch = Y_train[batch_index]
weights, bias, _, _ = opt.step(cost, weights, bias, feats_train_batch, Y_train_batch)
# Compute predictions on train and validation set
predictions_train = [np.sign(variational_classifier(weights, bias, f)) for f in feats_train]
predictions_val = [np.sign(variational_classifier(weights, bias, f)) for f in feats_val]
# Compute accuracy on train and validation set
acc_train = accuracy(Y_train, predictions_train)
acc_val = accuracy(Y_val, predictions_val)
print(
"Iter: {:5d} | Cost: {:0.7f} | Acc train: {:0.7f} | Acc validation: {:0.7f} "
"".format(it + 1, cost(weights, bias, features, Y), acc_train, acc_val)
)
Iter: 1 | Cost: 1.4647420 | Acc train: 0.4933333 | Acc validation: 0.5600000 Iter: 2 | Cost: 1.3817093 | Acc train: 0.4933333 | Acc validation: 0.5600000 Iter: 3 | Cost: 1.2571520 | Acc train: 0.4800000 | Acc validation: 0.5600000 Iter: 4 | Cost: 1.1050173 | Acc train: 0.4533333 | Acc validation: 0.5600000 Iter: 5 | Cost: 0.9704723 | Acc train: 0.4933333 | Acc validation: 0.6000000 Iter: 6 | Cost: 0.8841282 | Acc train: 0.5866667 | Acc validation: 0.7600000 Iter: 7 | Cost: 0.8459409 | Acc train: 0.7200000 | Acc validation: 0.7200000 Iter: 8 | Cost: 0.8408246 | Acc train: 0.7333333 | Acc validation: 0.6800000 Iter: 9 | Cost: 0.8453070 | Acc train: 0.6400000 | Acc validation: 0.6400000 Iter: 10 | Cost: 0.8410790 | Acc train: 0.6266667 | Acc validation: 0.6400000 Iter: 11 | Cost: 0.8081787 | Acc train: 0.6400000 | Acc validation: 0.6400000 Iter: 12 | Cost: 0.7677528 | Acc train: 0.7200000 | Acc validation: 0.6800000 Iter: 13 | Cost: 0.7381111 | Acc train: 0.8266667 | Acc validation: 0.7200000 Iter: 14 | Cost: 0.7194719 | Acc train: 0.8800000 | Acc validation: 1.0000000 Iter: 15 | Cost: 0.7206167 | Acc train: 0.8133333 | Acc validation: 0.9200000 Iter: 16 | Cost: 0.7321430 | Acc train: 0.7866667 | Acc validation: 0.8000000 Iter: 17 | Cost: 0.7440584 | Acc train: 0.7733333 | Acc validation: 0.7600000 Iter: 18 | Cost: 0.7489444 | Acc train: 0.7600000 | Acc validation: 0.7600000 Iter: 19 | Cost: 0.7472663 | Acc train: 0.7733333 | Acc validation: 0.7600000 Iter: 20 | Cost: 0.7414264 | Acc train: 0.7866667 | Acc validation: 0.7600000 Iter: 21 | Cost: 0.7340591 | Acc train: 0.8266667 | Acc validation: 0.9200000 Iter: 22 | Cost: 0.7235194 | Acc train: 0.8266667 | Acc validation: 0.9200000 Iter: 23 | Cost: 0.7159816 | Acc train: 0.8933333 | Acc validation: 0.9600000 Iter: 24 | Cost: 0.7110693 | Acc train: 0.8933333 | Acc validation: 1.0000000 Iter: 25 | Cost: 0.7064306 | Acc train: 0.8933333 | Acc validation: 1.0000000
plt.figure(figsize=(12,9))
cm = plt.cm.RdBu
# make data for decision regions
xx, yy = np.meshgrid(np.linspace(0.0, 1.5, 20), np.linspace(0.0, 1.5, 20))
X_grid = [np.array([x, y]) for x, y in zip(xx.flatten(), yy.flatten())]
# preprocess grid points like data inputs above
padding = 0.3 * np.ones((len(X_grid), 1))
X_grid = np.c_[np.c_[X_grid, padding], np.zeros((len(X_grid), 1))] # pad each input
normalization = np.sqrt(np.sum(X_grid ** 2, -1))
X_grid = (X_grid.T / normalization).T # normalize each input
features_grid = np.array(
[get_angles(x) for x in X_grid]
) # angles for state preparation are new features
predictions_grid = [variational_classifier(weights, bias, f) for f in features_grid]
Z = np.reshape(predictions_grid, xx.shape)
# plot decision regions
cnt = plt.contourf(
xx, yy, Z, levels=np.arange(-1, 1.1, 0.1), cmap=cm, alpha=0.8, extend="both"
)
plt.contour(
xx, yy, Z, levels=[0.0], colors=("black",), linestyles=("--",), linewidths=(0.8,)
)
plt.colorbar(cnt, ticks=[-1, 0, 1])
# plot data
plt.scatter(
X_train[:, 0][Y_train == 1],
X_train[:, 1][Y_train == 1],
c="b",
marker="o",
edgecolors="k",
label="class 1 train",
)
plt.scatter(
X_val[:, 0][Y_val == 1],
X_val[:, 1][Y_val == 1],
c="b",
marker="^",
edgecolors="k",
label="class 1 validation",
)
plt.scatter(
X_train[:, 0][Y_train == -1],
X_train[:, 1][Y_train == -1],
c="r",
marker="o",
edgecolors="k",
label="class -1 train",
)
plt.scatter(
X_val[:, 0][Y_val == -1],
X_val[:, 1][Y_val == -1],
c="r",
marker="^",
edgecolors="k",
label="class -1 validation",
)
plt.legend()
plt.show()
Here, we will represent simple model of Quantum convolutional neural network (QCNN), and compare it with classical convolutional neural network (CNN) model.
Classical convolutional neural network
If the input is an image, small local regions are sequentially processed with the same kernel. The results obtained for each region are usually associated to different channels of a single output pixel. The union of all the output pixels produces a new image-like object, which can be further processed by additional layers.
In summary, there are mainly three parts, convolutional layer, pooling, and classification.
convolutional layer: In a nutshell, convolution (합성곱) is a mathematical operator where taking many data values and combine them into different values. "Highlighting the valuable part of the information."pooling: "Extracting valuable information from the data, and then resizing them into smaller size."Classification: At the end, we apply the neural network model (신경망 모형) to classify 'filtered' data into given classes.Quantum convolutional neural network
Just like classical convolutional neural network, we can apply the similar architecture of the model into quantum circuit. A simple scheme is also represented in the figure at the top.
The main difference with respect to a classical convolution is that a quantum circuit can generate highly complex kernels whose computation could be, at least in principle, classically intractable.